# Beta gradient assignment

#### Michal Kowalewski

#### August 2, 2018

Beta gradient analysis is explictly spatial. We will examine a dataset of marine invertebrate samples from coastal North Carolina. First, read and examine live community data and associated locality info.

We will also need the library ‘vegan’.

    library(vegan)

    ## Loading required package: permute

    ## Loading required package: lattice

    ## This is vegan 2.5-2

    live <- read.csv('Final_All_Taxa_BetaDiv.csv',head=TRUE)
    refs <- read.csv('Refs_local.csv')

Now, assess beta diversity gradient as a function of depth difference between localities. Bray-Curtis is an obvious (safe) measure, but you can use other measures as well.

The obvious way (although perhaps not CPU-efficient) is to use nested for loops. Try those first.

    # for loop solution
    out1 <- NULL
    for (i in 1:nrow(live)) {
     for (j in 1:nrow(live)) {
      if (j > i) {
      dist <- 1 - vegdist(rbind(live[i,], live[j,]))
      D.depth <- abs(refs$Depth[i] - refs$Depth[j])
      habitat <- ifelse(refs$Habitat[i] == refs$Habitat[j],
                        habitat <- 1, habitat <- 0)
      out1 <- rbind(out1, c(dist=dist, depth=D.depth, hab=habitat,
                          min.n=min(sum(live[i,]), sum(live[j,]))))
      }
     }
    }
    out3 <- data.frame(out1)

You may try to utilize a well tested and efficient function such as vegdist to see if you can make your script faster.

    # making use of vegdist function
    min.nveg <- cbind(designdist(live, method = "A", terms='minimum'),
                      designdist(live, method = "B", terms='minimum'))
    min.nveg2 <- apply(min.nveg, 1, min)
    habdist <- vegdist(refs$Habitat, 'euclidean')
    hab1 <- ifelse(habdist==0, hab1 <- 1, hab1 <- 0)
    out2 <- cbind(dist = 1 - vegdist(live),
                  depth = vegdist(refs$Depth, 'euclidean'),
                  hab = hab1, min.n=min.nveg2)
    out4 <- data.frame(out2)

If you created multiple ways to do this exercise, check if those methods agree. Say, your outputs are ‘out3’ and ‘out4’. Check if they are the same.

    sum(out4!=out3) # check if all values agree

    ## [1] 0

The final step is to visualize your results. If you generated multiple solutions. It may be best to write a custom functioon for plotting your data.

    bg.plot <- function(out) {
    bubble <- (out$min.n - min(out$min.n))/(max(out$min.n) - min(out$min.n)) + 1
      b1 <- tapply(out$dist, round(out$depth, -0.5), mean)
      b2 <- tapply(out$dist[out$hab==1], round(out$depth[out$hab==1], -0.5), mean)
      b3 <- tapply(out$dist[out$hab==0], round(out$depth[out$hab==0], -0.5), mean)
    plot(out[,2], out[,1], ylim=c(0,1), pch=21, cex=bubble, las=1,
         col=c('tan3', 'skyblue')[as.factor(out$hab)],
         bg=adjustcolor(c('tan3', 'skyblue')[as.factor(out$hab)], 0.3),
         xlab=bquote(Delta~'depth [m]'),
         ylab='Bray-Curtis similarity')
       points(as.numeric(names(b1)), b1, type='l', col='white', lwd=5)
       points(as.numeric(names(b1)), b1, type='l', col='black', lwd=2.5)
       points(as.numeric(names(b2)), b2, type='l', col='white', lwd=5)
       points(as.numeric(names(b2)), b2, type='l', col='tan3', lwd=2.5)
       points(as.numeric(names(b3)), b3, type='l', col='white', lwd=5)
       points(as.numeric(names(b3)), b3, type='l', col='skyblue', lwd=2.5)
      legend('topright', lty=1, col=c('black', 'skyblue', 'tan3'),
           c('all', 'between habitats', 'within-habitat'), cex=0.7)
      legend('topleft', pch=21, pt.cex=c(max(bubble), min(bubble)),
           legend=c(paste('largest n =', max(out$min.n)),
                    paste('smallest n =', min(out$min.n))))
    }

Now, plot your outcomes.

#### Solution 1

    bg.plot(out3)

#### Solution 2

    bg.plot(out4)

#### adding turnover model (exponential fit)

This function has not been tested(use at your own peril!)

    decay<-function(x){
      fitH<- x[which(x[,1]*x[,2]>0),] # get rid of 0 values (not allowed for log operations)
      expH<- lm(log(fitH[,2])~ fitH[,1]) # estimate exponent
      xlH<- seq(0.1,15,0.1) # set a series of points for range of depth difference values (for plotting over a given range of x values)
      ylH<- exp(expH$coefficients[2]*xlH+expH$coefficients[1]) # compute predicted y values (bray-curtis values)
      adjr2 <- summary(expH)$adj.r.squared # compute adjusted r-square
      decay<-cbind(xlH,ylH,expH$coefficients[2])
      return(decay)
    }

Now plot the plot again but add a decay model line

    bg.plot(out4)
     points(decay(cbind(out4$depth, out4$dist)), type='l', col='white', lwd=4)
     points(decay(cbind(out4$depth, out4$dist)), type='l', col='red1', lwd=2)

And then cite packages…

    citation('vegan')

    ##
    ## To cite package 'vegan' in publications use:
    ##
    ##   Jari Oksanen, F. Guillaume Blanchet, Michael Friendly, Roeland
    ##   Kindt, Pierre Legendre, Dan McGlinn, Peter R. Minchin, R. B.
    ##   O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens,
    ##   Eduard Szoecs and Helene Wagner (2018). vegan: Community Ecology
    ##   Package. R package version 2.5-2.
    ##   https://CRAN.R-project.org/package=vegan
    ##
    ## A BibTeX entry for LaTeX users is
    ##
    ##   @Manual{,
    ##     title = {vegan: Community Ecology Package},
    ##     author = {Jari Oksanen and F. Guillaume Blanchet and Michael Friendly and Roeland Kindt and Pierre Legendre and Dan McGlinn and Peter R. Minchin and R. B. O'Hara and Gavin L. Simpson and Peter Solymos and M. Henry H. Stevens and Eduard Szoecs and Helene Wagner},
    ##     year = {2018},
    ##     note = {R package version 2.5-2},
    ##     url = {https://CRAN.R-project.org/package=vegan},
    ##   }
    ##
    ## ATTENTION: This citation information has been auto-generated from
    ## the package DESCRIPTION file and may need manual editing, see
    ## 'help("citation")'.